******************************************************************************
*										
*                           Boese & Eberhardt							
*            
*                          Facets of Democracy
*													
******************************************************************************
*                   PCDID - heterogeneous model (1949-2018)
*			    (Specification with exports/trade & pop growth)
* 				  Low-level indices: Interaction Effects
******************************************************************************
*                       Using standardized thresholds 
******************************************************************************
*                           Created: 17th July 2024					
*                      Last Changed: 18th July 2024					
******************************************************************************
*                              Markus Eberhardt						
******************************************************************************
*                 School of Economics, University of Nottingham,			
*                  University Park, Nottingham NG7 2RD, England			
*                     email: markus.eberhardt@nottingham.ac.uk				
*                   web: https://sites.google.com/site/medevecon/			
******************************************************************************

clear all
set matsize 11000

***
*** Paths
***

global path "D:/Dropbox/Facets of Democracy"
global rawpath "D:/Dropbox"
global path "/Users/lezme/Dropbox/Facets of Democracy"
global rawpath "/Users/lezme/Dropbox"

global data "$path/Data"
global dofolder "$path/Do_files"
global otherdata "$path/Data"
global outputfolder "$path/Output"
global texfolder "$rawpath/Apps/Overleaf/Facets of Democracy"

global xtmgpath = "D:/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"
global xtmgpath = "/Users/lezme/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"

****
**** Get merged dataset and prepare variables
****

use "$data/madd_ert_vdem_2021.dta", clear
tsset nwbcode year

* Maddison: per capita GDP (in logs)*100 and population growth rate
gen y = ln(gdppc)*100
gen lpop = ln(pop)
gen dlpop = 100*d.lpop

* DOTS data: export/trade
gen ex = ex_trade*100

xtreg y dlpop ex v2x_libdem if year>=1959, fe

drop csum 
gen c=1 if e(sample)
by nwbcode: egen csum=sum(c) if c==1
gen sample= (e(sample))

* V-Dem: democracy dummies from high-, mid- and low-level indicators
local z "_liberal _polyarchy _freexp_altinf _frassoc_thick el_frefair cl_rol _jucon lg_legcon"
foreach k of local z{
	sum v2x`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st`k' = (v2x`k'-mean`k')/sd`k' if sample==1
}


* Positive values
local z "0 25 50"
foreach k of local z{
	local norm_`k' = `k'/100
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal)
	gen dem_exp_`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legcon_`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legcon_`k' "Leg. Constraints Threshold > mean + .`k' sd"
}

local z "125"
foreach k of local z{
	local norm_`k' = `k'/1000
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal)
	gen dem_exp_`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legcon_`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legcon_`k' "Leg. Constraints Threshold > mean + .`k' sd"
}

* Negative values
local z "25 50"
foreach k of local z{
	local norm_`k' = -`k'/100
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	gen dem_lib_neg`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal)
	gen dem_exp_neg`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_neg`k' "Freedom of Expression Threshold > mean - .`k' sd"
	gen dem_assoc_neg`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_neg`k' "Freedom of Association Threshold > mean - .`k' sd"
	gen dem_fairel_neg`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean - .`k' sd"
	gen dem_rolaw_neg`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_neg`k' "Rule of Law Threshold > mean - .`k' sd"
	gen dem_jucon_neg`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_neg`k' "Jud. Constraints Threshold > mean - .`k' sd"
	gen dem_legcon_neg`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legcon_neg`k' "Leg. Constraints Threshold > mean - .`k' sd"
}

local z "125"
foreach k of local z{
	local norm_`k' = -`k'/1000
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	gen dem_lib_neg`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal)
	gen dem_exp_neg`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_neg`k' "Freedom of Expression Threshold > mean - .`k' sd"
	gen dem_assoc_neg`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_neg`k' "Freedom of Association Threshold > mean - .`k' sd"
	gen dem_fairel_neg`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean - .`k' sd"
	gen dem_rolaw_neg`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_neg`k' "Rule of Law Threshold > mean - .`k' sd"
	gen dem_jucon_neg`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_neg`k' "Jud. Constraints Threshold > mean - .`k' sd"
	gen dem_legcon_neg`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legcon_neg`k' "Leg. Constraints Threshold > mean - .`k' sd"
}

save "$data/madd_ert_vdem_2021_helper_inter_low.dta", replace
 
****
**** PCDID 
**** 

estimates clear 
cls 

local i = 1

**** Select the variables, sample period, factor-augmentation 
global xvars1 		"ex dlpop" // X variables used to construct factors 
global xvars2 		"ex dlpop" // Xs included in the PCDID (leave blank for 'plain vanilla' model)
global factors 		"6"		 // max # of factors (over-prediction is not as bad as underprediction)
global uhat1 		"uhatMG1" // uhatMG: residuals from heterog regression (first control group)
global uhat2 		"uhatMG2" // uhatMG: residuals from heterog regression (second control group)
global uhatint 		"uhatMGint" // uhatMG: residuals from heterog regression (interaction control group)

*** Pick which threshold definition is used (mean or mean +- sd)
*local sdvari "neg25 neg125 0 125 25"
local sdvari "neg25 neg125 0 125 25"
* 0 is the model using the standardised mean of the democracy index
* the other indices are neg25 neg125 125 25
foreach ksdvari of local sdvari{
		
	*** List democracy indicators employed here (list1 always poly or lib; list_1 can be several lower-level indicators)
	
	*global demo_list1	"poly_`ksdvari'" 
	global demo_list1	"lib_`ksdvari'" 
	*global demo_list2	"dem_poly_`ksdvari'" 
	global demo_list2	"dem_lib_`ksdvari'" 

	*global demo_list_1	"rolaw_`ksdvari' jucon_`ksdvari' legcon_`ksdvari'" 
	global demo_list_1	"exp_`ksdvari' assoc_`ksdvari' fairel_`ksdvari'" 
	*global demo_list_2	"dem_exp_dem_rolaw_`ksdvari' dem_jucon_`ksdvari' dem_legcon_`ksdvari'" 
	global demo_list_2	"dem_exp_`ksdvari' dem_assoc_`ksdvari' dem_fairel_`ksdvari'" 

 	run "$xtmgpath/xtmg.ado"

	* Start of loop for different end years 
	local zzz "2018"
	foreach kkk of local zzz{
		global endyear 		"`kkk'"		// 2018

		* Start of loop for different start years 
		local zz "1959"
		foreach kk of local zz{
			global startyear 	"`kk'"	// 1949

			* Start of loop for the first democracy indicators (e.g. poly) 
			local z "$demo_list1"
			foreach k of local z{
				global democracy 	"`k'"	// picks one democracy indicator

				* Start of loop for the second democracy indicator (e.g. rolaw) !!! Underscore_1 !!!
				local z2 "$demo_list_1"
				foreach k2 of local z2{
					global democracy_too 	"`k2'"	// picks second democracy indicator (from a list)

					* Start of factor loop: how many factors for dem_list1 variable to be included in the PCDID
					forvalues f=1/$factors{
					preserve
					quietly{
						
						tsset nwbcode year
						keep if year>=$startyear & year<=$endyear
						
						gen dem_dum = dem_$democracy
						gen dem_dumm = dem_$democracy_too
						gen dem_inter = dem_dum * dem_dumm
						sort nwbcode year
						
						* Cross the threshold
						gen count_$democracy = 1 if ((dem_dum==0 & l.dem_dum==1) | (dem_dum==1 & l.dem_dum==0))
						gen count_$democracy_too = 1 if ((dem_dumm==0 & l.dem_dumm==1) | (dem_dumm==1 & l.dem_dumm==0))
						gen count_inter = count_$democracy * count_$democracy_too
						
						by nwbcode: egen dem_dum_count_$democracy = sum(count_$democracy)
						by nwbcode: egen dem_dumm_count_$democracy_too = sum(count_$democracy_too)
						by nwbcode: egen dem_dum_count_inter = sum(count_inter)
						
						replace count_$democracy = 0 if (dem_dum==0 & l.dem_dum==1)
						replace count_$democracy_too = 0 if (dem_dumm==0 & l.dem_dumm==1)
						replace count_inter = 0 if ((dem_dumm==0 & l.dem_dumm==1) | dem_dum==0 & l.dem_dum==1)
						
						* Cross the threshold from below
						by nwbcode: egen trudem_dum_count_$democracy = sum(count_$democracy)
						by nwbcode: egen trudem_dumm_count_$democracy_too = sum(count_$democracy_too)
						by nwbcode: egen trudem_dum_count_inter = sum(count_inter)

						by nwbcode: egen dem_dum_exposure_$democracy = sum(dem_dum)
						by nwbcode: egen dem_dumm_exposure_$democracy_too = sum(dem_dumm)
						by nwbcode: egen dem_dum_exposure_inter = sum(dem_inter)
	
						by nwbcode: egen sample_obs_$democracy = sum(sample)
						by nwbcode: egen sample_obs_$democracy_too = sum(sample)
						by nwbcode: egen sample_obs_inter = sum(sample)

						gen autoc_exposure_$democracy = sample_obs_$democracy - dem_dum_exposure_$democracy
						gen autoc_exposure_$democracy_too = sample_obs_$democracy_too - dem_dumm_exposure_$democracy_too
						gen autoc_exposure_inter = sample_obs_inter - dem_dum_exposure_inter

						*** 1 ***
						* Sample identifiers for the first dem variable
						xtlogit dem_dum y if sample==1, fe iterate(1)
						* Treat are those which experienced crossing the threshold
						gen treat = (e(sample))
						replace treat = . if sample==0
						* Below are those which are always below the threshold
						gen below = (!e(sample) & sample==1 & !missing(dem_$democracy))
						replace below = . if sample==0
						* Always above threshold (adjust)
						gen above = (below==1 & dem_$democracy ==1)
						replace above = . if sample==0
						* Always below threshold (adjust)
						replace below = . if below==1 & dem_$democracy==1
						
						*** 2 ***
						* Sample identifiers for the second dem variable
						xtlogit dem_dumm y if sample==1, fe iterate(1)
						* Treat are those which experienced crossing the threshold
						gen treat_too = (e(sample))
						replace treat_too = . if sample==0
						* Below are those which are always below the threshold
						gen below_too = (!e(sample) & sample==1 & !missing(dem_$democracy_too))
						replace below_too = . if sample==0
						* Always above threshold (adjust)
						gen above_too = (below_too==1 & dem_$democracy_too ==1)
						replace above_too = . if sample==0
						* Always below threshold (adjust)
						replace below_too = . if below_too==1 & dem_$democracy_too==1

						*** 3 ***
						* Sample identifiers for the interaction variable
						gen treat_inter = treat_too*treat
						gen below_inter = below_too*below
						gen above_inter = above_too*above
						
						/*
						EXCLUDE: ANALYSIS OF THE TWO ELEMENTS OF THE INTERACTION
						*** 1 ***
						*** Estimate linear heterog regression of y on x AND second dem dummy in the first control group and collect residuals
						xtmg y dem_dumm $xvars1 if below==1 & sample==1 & above==0, res(uhatMG1)
						local below = e(N_g)
						local belowobs = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y dem_dumm $xvars1 if below==1 & sample==1 & above==0, fe
						predict uhatP1, e
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhat1 if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhat1 if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhat1 if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c1factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 1: c1factor*N
						
						*** 2 ***
						*** Estimate linear heterog regression of y on x AND first dem dummy in the second control group and collect residuals
						xtmg y dem_dum $xvars1 if below_too==1 & sample==1 & above_too==0, res(uhatMG2)
						local below_too = e(N_g)
						local belowobs_too = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y dem_dum $xvars1 if below_too==1 & sample==1  & above_too==0, fe
						predict uhatP2, e
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhat2 if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhat2 if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhat2 if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c2factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 2: c2factor*N

						*/
						
						*** 3 ***
						*** Estimate linear heterog regression of y on x in the interaction control group and collect residuals
						xtmg y $xvars1 if below_inter==1 & sample==1 & above_inter==0, res(uhatMGint)
						local below_inter = e(N_g)
						local belowobs_inter = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y $xvars1 if below_inter==1 & sample==1  & above_inter==0, fe
						predict uhatPint, e
						*** Compute residual period average in control group 
						sort year 
						by year: egen uCT = mean($uhatint)
						sort nwbcode year					
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhatint if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhatint if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhatint if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c3factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 3 for interaction: c3factor*N
						
						gen dem_interaction = dem_$democracy * dem_$democracy_too
						
						*** Test of controls ***
						xtlogit dem_interaction y $xvars2 ///
							if sample==1 & treat==1, fe iterate(1)
						xtmg dem_interaction $xvars2 c3factor*N ///
							if e(sample)
						testparm $xvars2
						local chi2_stat = r(chi2)
						local chi2_p = r(p)

						*** Alpha test auxiliary regressions ***
						xtmg y uCT ///
							dem_interaction ///
							$xvars2 ///
							if sample==1 & treat_inter==1
						mat alphas = e(betas)
						svmat alphas, name(alphas)
						gen alpha01 = alphas1-1
						replace alpha01 = . if (alphas2)==0
						mkmat alpha01, mat(alpha01)
						set seed 280115
						reg alpha01
						mat a1 = e(b)
						mat vv1 = e(V)
						scalar alphaest1 = a1[1,1]
						local alpha_est1 = alphaest1
						scalar vv1s = vv1[1,1]
						scalar se1s = sqrt(vv1s)
						scalar alphat1 = alphaest1/se1s
						local alpha_tstat1 = alphat1

						*** Run PCDID model ***
						eststo: xtmg y dem_interaction $xvars2 ///
							c3factor*N ///
							if sample==1 & treat_inter==1, res(resid) robust
						mat betas = e(betas)
						estadd scalar treated = e(N_g)
						estadd scalar obs = e(N)
						local treated = e(N_g)
						local obs = e(N)
						mat b1 = e(b)
						mat v1 = e(V)
						scalar b_1 = b1[1,1]
						scalar se_1 = sqrt(v1[1,1])	
						scalar mse = e(sigma)
						estadd scalar mse = e(sigma)
						capture xtcd resid if e(sample)
						capture scalar cdtest = r(pesaran)
						capture estadd scalar cdtest = r(pesaran)
						estadd scalar fac = `f'
						estadd scalar x_test_chi2 = `chi2_stat'
						estadd scalar x_test_p = `chi2_p'
						estadd scalar alpha_b1 = `alpha_est1'
						estadd scalar alpha_t1 = `alpha_tstat1'
						*estadd scalar controlN1 = `below'
						*estadd scalar controlobs1 = `belowobs'
						*estadd scalar controlN2 = `below_too'
						*estadd scalar controlobs2 = `belowobs_too'
						estadd scalar controlN3 = `below_inter'
						estadd scalar controlobs3 = `belowobs_inter'
						estadd scalar final_year = `kkk'
						estadd scalar start_year = `kk'
						est store model`i'
						
						keep if sample==1 & treat_inter==1
						tabstat nwbcode if sample==1 & treat_inter==1, by(wbcode) save
						tabstatmat nwbcodestat, nototal
						tabstat year if sample==1 & treat_inter==1, by(wbcode) stats(min max) save
						tabstatmat yearstat, nototal
						tabstat dem_dum_count_inter trudem_dum_count_inter dem_dum_exposure_inter ///
							sample_obs_inter autoc_exposure_inter, by(wbcode) stats(mean) save
						tabstatmat exposure, nototal
						drop _all
						
						svmat nwbcodestat, name(nwbcode)
						rename nwbcode1 nwbcode
						svmat yearstat, name(year)
						rename year1 startyear
						rename year2 endyear
						svmat exposure, name(count)
						rename count1 flips
						rename count2 democratisations
						rename count3 exposure
						rename count4 total_obs
						rename count5 nonexposure
						svmat betas, name(beta)
						gen wbcode=""
						replace wbcode="AFG" if nwbcode==	2
						replace wbcode="AGO" if nwbcode==	3
						replace wbcode="ALB" if nwbcode==	5
						replace wbcode="ARE" if nwbcode==	7
						replace wbcode="ARG" if nwbcode==	8
						replace wbcode="ARM" if nwbcode==	9
						replace wbcode="AUS" if nwbcode==	12
						replace wbcode="AUT" if nwbcode==	13
						replace wbcode="AZE" if nwbcode==	14
						replace wbcode="BDI" if nwbcode==	15
						replace wbcode="BEL" if nwbcode==	17
						replace wbcode="BEN" if nwbcode==	18
						replace wbcode="BFA" if nwbcode==	19
						replace wbcode="BGD" if nwbcode==	20
						replace wbcode="BGR" if nwbcode==	21
						replace wbcode="BHR" if nwbcode==	22
						replace wbcode="BIH" if nwbcode==	24
						replace wbcode="BLR" if nwbcode==	25
						replace wbcode="BOL" if nwbcode==	29
						replace wbcode="BRA" if nwbcode==	30
						replace wbcode="BRB" if nwbcode==	31
						replace wbcode="BWA" if nwbcode==	36
						replace wbcode="CAF" if nwbcode==	37
						replace wbcode="CAN" if nwbcode==	38
						replace wbcode="CHE" if nwbcode==	39
						replace wbcode="CHL" if nwbcode==	40
						replace wbcode="CHN" if nwbcode==	41
						replace wbcode="CIV" if nwbcode==	42
						replace wbcode="CMR" if nwbcode==	43
						replace wbcode="COG" if nwbcode==	45
						replace wbcode="COL" if nwbcode==	46
						replace wbcode="COM" if nwbcode==	47
						replace wbcode="CPV" if nwbcode==	48
						replace wbcode="CRI" if nwbcode==	49
						replace wbcode="CUB" if nwbcode==	51
						replace wbcode="CYP" if nwbcode==	52
						replace wbcode="CZE" if nwbcode==	53
						replace wbcode="DEU" if nwbcode==	55
						replace wbcode="DJI" if nwbcode==	56
						replace wbcode="DNK" if nwbcode==	58
						replace wbcode="DOM" if nwbcode==	59
						replace wbcode="DZA" if nwbcode==	60
						replace wbcode="ECU" if nwbcode==	61
						replace wbcode="EGY" if nwbcode==	62
						replace wbcode="ESP" if nwbcode==	64
						replace wbcode="EST" if nwbcode==	65
						replace wbcode="ETH" if nwbcode==	66
						replace wbcode="FIN" if nwbcode==	67
						replace wbcode="FRA" if nwbcode==	69
						replace wbcode="GAB" if nwbcode==	71
						replace wbcode="GBR" if nwbcode==	72
						replace wbcode="GEO" if nwbcode==	73
						replace wbcode="GHA" if nwbcode==	74
						replace wbcode="GIN" if nwbcode==	76
						replace wbcode="GMB" if nwbcode==	77
						replace wbcode="GNB" if nwbcode==	78
						replace wbcode="GNQ" if nwbcode==	79
						replace wbcode="GRC" if nwbcode==	80
						replace wbcode="GTM" if nwbcode==	83
						replace wbcode="HKG" if nwbcode==	87
						replace wbcode="HND" if nwbcode==	89
						replace wbcode="HRV" if nwbcode==	91
						replace wbcode="HTI" if nwbcode==	92
						replace wbcode="HUN" if nwbcode==	93
						replace wbcode="IDN" if nwbcode==	95
						replace wbcode="IND" if nwbcode==	96
						replace wbcode="IRL" if nwbcode==	97
						replace wbcode="IRN" if nwbcode==	98
						replace wbcode="IRQ" if nwbcode==	99
						replace wbcode="ISL" if nwbcode==	100
						replace wbcode="ISR" if nwbcode==	101
						replace wbcode="ITA" if nwbcode==	102
						replace wbcode="JAM" if nwbcode==	103
						replace wbcode="JOR" if nwbcode==	104
						replace wbcode="JPN" if nwbcode==	105
						replace wbcode="KAZ" if nwbcode==	106
						replace wbcode="KEN" if nwbcode==	107
						replace wbcode="KGZ" if nwbcode==	108
						replace wbcode="KHM" if nwbcode==	109
						replace wbcode="KOR" if nwbcode==	112
						replace wbcode="KWT" if nwbcode==	113
						replace wbcode="LAO" if nwbcode==	115
						replace wbcode="LBN" if nwbcode==	116
						replace wbcode="LBR" if nwbcode==	117
						replace wbcode="LBY" if nwbcode==	118
						replace wbcode="LKA" if nwbcode==	120
						replace wbcode="LSO" if nwbcode==	121
						replace wbcode="LTU" if nwbcode==	122
						replace wbcode="LUX" if nwbcode==	123
						replace wbcode="LVA" if nwbcode==	124
						replace wbcode="MAR" if nwbcode==	126
						replace wbcode="MDA" if nwbcode==	128
						replace wbcode="MDG" if nwbcode==	129
						replace wbcode="MEX" if nwbcode==	132
						replace wbcode="MLI" if nwbcode==	135
						replace wbcode="MLT" if nwbcode==	136
						replace wbcode="MMR" if nwbcode==	137
						replace wbcode="MNE" if nwbcode==	138
						replace wbcode="MNG" if nwbcode==	139
						replace wbcode="MOZ" if nwbcode==	140
						replace wbcode="MRT" if nwbcode==	141
						replace wbcode="MUS" if nwbcode==	143
						replace wbcode="MWI" if nwbcode==	144
						replace wbcode="MYS" if nwbcode==	145
						replace wbcode="NAM" if nwbcode==	146
						replace wbcode="NER" if nwbcode==	148
						replace wbcode="NGA" if nwbcode==	149
						replace wbcode="NIC" if nwbcode==	150
						replace wbcode="NLD" if nwbcode==	151
						replace wbcode="NOR" if nwbcode==	152
						replace wbcode="NPL" if nwbcode==	153
						replace wbcode="NZL" if nwbcode==	156
						replace wbcode="OMN" if nwbcode==	158
						replace wbcode="PAK" if nwbcode==	159
						replace wbcode="PAN" if nwbcode==	160
						replace wbcode="PER" if nwbcode==	161
						replace wbcode="PHL" if nwbcode==	162
						replace wbcode="POL" if nwbcode==	165
						replace wbcode="PRK" if nwbcode==	168
						replace wbcode="PRT" if nwbcode==	170
						replace wbcode="PRY" if nwbcode==	171
						replace wbcode="QAT" if nwbcode==	176
						replace wbcode="RUS" if nwbcode==	179
						replace wbcode="RWA" if nwbcode==	180
						replace wbcode="SAU" if nwbcode==	181
						replace wbcode="SDN" if nwbcode==	183
						replace wbcode="SEN" if nwbcode==	184
						replace wbcode="SGP" if nwbcode==	185
						replace wbcode="SLE" if nwbcode==	187
						replace wbcode="SLV" if nwbcode==	188
						replace wbcode="STP" if nwbcode==	195
						replace wbcode="SVK" if nwbcode==	198
						replace wbcode="SVN" if nwbcode==	199
						replace wbcode="SWE" if nwbcode==	200
						replace wbcode="SWZ" if nwbcode==	201
						replace wbcode="SYC" if nwbcode==	203
						replace wbcode="SYR" if nwbcode==	204
						replace wbcode="TCD" if nwbcode==	205
						replace wbcode="TGO" if nwbcode==	206
						replace wbcode="THA" if nwbcode==	207
						replace wbcode="TJK" if nwbcode==	208
						replace wbcode="TKM" if nwbcode==	209
						replace wbcode="TTO" if nwbcode==	214
						replace wbcode="TUN" if nwbcode==	215
						replace wbcode="TUR" if nwbcode==	216
						replace wbcode="TZA" if nwbcode==	220
						replace wbcode="UGA" if nwbcode==	221
						replace wbcode="UKR" if nwbcode==	222
						replace wbcode="URY" if nwbcode==	223
						replace wbcode="USA" if nwbcode==	224
						replace wbcode="UZB" if nwbcode==	225
						replace wbcode="VEN" if nwbcode==	228
						replace wbcode="VNM" if nwbcode==	229
						replace wbcode="YEM" if nwbcode==	235
						replace wbcode="ZAF" if nwbcode==	238
						replace wbcode="ZMB" if nwbcode==	240
						replace wbcode="ZWE" if nwbcode==	241
						order nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
						keep nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
						gen regime_def = "inter"
						local startyearx "$startyear"
						local endyearx "$endyear" 
						local democracyx "$democracy" 
						local democracyxx "$democracy_too" 
						save "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_`startyearx'_`endyearx'_`democracyx'_`democracyxx'_`f'.dta", replace
						
						local i = `i'+1

					}
					* end of quietly
												
					di in ye _n "*****************************************************************"
					di in ye "*** Single regime cutoff defined by " in gr "dem_$democracy" 
					di in ye "*****************************************************************"
					di in ye "*** PCDID (MG) estimate with" in gr " `f' " in ye "common factor(s) included      ***"
					di in ye "*****************************************************************"
					di in ye "Nominal sample period: " in gr "$startyear" in ye " to " in gr "$endyear "
					di in ye "resulting in " in gr "`treated'" in ye " treated countries with " ///
						in gr "`obs'" in ye " observations. "
					di in ye "MSE: " in gr %5.3f mse
					*di in ye "Control sample 1: " in gr "`below'" in ye " countries with " in gr "`belowobs'" in ye " obs."
					*di in ye "Control sample 2: " in gr "`below_too'" in ye " countries with " in gr "`belowobs_too'" in ye " obs."
					di in ye "Control sample: " in gr "`below_inter'" in ye " countries with " in gr "`belowobs_inter'" in ye " obs."
					di in ye "*****************************************************************"
					*di in ye "Estimates for 1st democracy indicator: " in gr "dem_$democracy"
					*di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
					*	in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
					*di in ye "*****************************************************************"
					di in ye "Estimates for interaction with: "  in gr "dem_$democracy_too"
					di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
						in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
					di in ye "*****************************************************************"
					di in ye "*** Alpha Test: " in gr  %5.3f alphaest1 in ye " t-stat: " in gr  %5.3f alphat1
					di in ye "*****************************************************************"
					di in ye "*** Test of controls " in gr %5.2f `chi2_stat' in ye " (" in gr %5.3f `chi2_p' in ye ").      ***"
					di in ye "*****************************************************************"
					local ii = `i'-1
					di in red "`ii'"
					
					restore	
					
					}
					* end of the factor loop
				}
				* end of second dem indicator loop
			}
			* end of first dem indicator loop 
		}
		* end of start year loop
	}
	* end of end year loop

	local mytitle "PCDID"
	local esttabformat  "nodepvars lines compress label nogaps numbers star(* 0.10 ** 0.05 *** 0.01)"
	local esttabstats "stats(start_year final_year fac treated obs controlN1 controlobs1 controlN2 controlobs2 controlN3 controlobs3 x_test_chi2 x_test_p alpha_b1 alpha_t1, fmt(0 0 0 0 0 0 0 0 0 2 2 3 3))"
	local esttabkeep "keep($demo_list2 $demo_list_2 dem_interaction $xvars2) order($demo_list2 $demo_list_2 dem_interaction $xvars2)"
	esttab model* using "$outputfolder/20240717_Maddison2021_PCDID_Drilling_`democracyx'_interaction_simple_all_factors.csv", b(3) se(3) abs br ///
		`esttabformat' `esttabkeep' `esttabstats' style(tab) replace nodep

}
* end of the loop related to the variation of the threshold (-.25sd -.125sd 0 [mean] +.125sd +.25sd)


exit

***
*** Graph: Comparison interaction with poly
***

use "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_poly_0_diagnostics_4.dta", clear 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_poly
}
drop nwbcode
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_0_rolaw_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_polyrolaw
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_0_jucon_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_polyjucon
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_0_legcon_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_polylegcon
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_rolaw_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_rolaw
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_jucon_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_jucon
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_legon_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_legcon
}
drop nwbcode _merge

*** No horserace

*** POLY & ROLAW
tab startyear_poly, gen(by_)
mrunning beta_poly exposure_poly flips_poly if !missing(beta_poly,beta_rolaw, beta_polyrolaw), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_poly
rename changese__1 changese_poly
drop change__* changese__*
gen tstat_poly = change_poly/changese_poly
drop by_*

tab startyear_rolaw, gen(by_)
mrunning beta_rolaw exposure_rolaw flips_rolaw if !missing(beta_poly,beta_rolaw, beta_polyrolaw), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_rolaw
rename changese__1 changese_rolaw
drop change__* changese__*
gen tstat_rolaw = change_rolaw/changese_rolaw
drop by_*

tab startyear_polyrolaw, gen(by_)
mrunning beta_polyrolaw exposure_polyrolaw flips_polyrolaw if !missing(beta_poly,beta_rolaw, beta_polyrolaw), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_polyrolaw
rename changese__1 changese_polyrolaw
drop change__* changese__*
gen tstat_polyrolaw = change_polyrolaw/changese_polyrolaw
drop by_*

* 68 countries 

twoway ///
	(hist exposure_polyrolaw, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_poly exposure_poly, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(line change_rolaw exposure_rolaw, sort ///
		lw(.8) lcolor(gs7*1) lpattern(dash)) ///
	(scatter change_rolaw exposure_rolaw if abs(tstat_rolaw)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_rolaw exposure_rolaw if abs(tstat_rolaw)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(line change_polyrolaw exposure_polyrolaw, sort ///
		lw(.8) lcolor(midblue*1) lpattern(solid)) ///
	(scatter change_polyrolaw exposure_polyrolaw if abs(tstat_polyrolaw)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(midblue*1.25%60) mlwidth(thin)) ///
	(scatter change_polyrolaw exposure_polyrolaw if abs(tstat_polyrolaw)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(midblue%60) mlcolor(midblue*1.25%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-20.1 40.1))   ///
	ylabel(-20(10)40, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)12, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "68") label(2 "Polyarchy > mean") label(5 "Rule of Law > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_poly_rolaw.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_poly_rolaw.eps", as(eps) replace

drop change_poly changese_poly tstat_poly

*** POLY & JUCON
tab startyear_poly, gen(by_)
mrunning beta_poly exposure_poly flips_poly if !missing(beta_poly,beta_jucon, beta_polyjucon), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_poly
rename changese__1 changese_poly
drop change__* changese__*
gen tstat_poly = change_poly/changese_poly
drop by_*

tab startyear_jucon, gen(by_)
mrunning beta_jucon exposure_jucon flips_jucon if !missing(beta_poly,beta_jucon, beta_polyjucon), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_jucon
rename changese__1 changese_jucon
drop change__* changese__*
gen tstat_jucon = change_jucon/changese_jucon
drop by_*

tab startyear_polyjucon, gen(by_)
mrunning beta_polyjucon exposure_polyjucon flips_polyjucon if !missing(beta_poly,beta_jucon, beta_polyjucon), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_polyjucon
rename changese__1 changese_polyjucon
drop change__* changese__*
gen tstat_polyjucon = change_polyjucon/changese_polyjucon
drop by_*
* 51 countries 

twoway ///
	(hist exposure_polyjucon, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_poly exposure_poly, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(line change_jucon exposure_jucon, sort ///
		lw(.8) lcolor(gs7*1) lpattern(dash)) ///
	(scatter change_jucon exposure_jucon if abs(tstat_jucon)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_jucon exposure_jucon if abs(tstat_jucon)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(line change_polyjucon exposure_polyjucon, sort ///
		lw(.8) lcolor(pink*1) lpattern(solid)) ///
	(scatter change_polyjucon exposure_polyjucon if abs(tstat_polyjucon)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(pink*1.5%60) mlwidth(thin)) ///
	(scatter change_polyjucon exposure_polyjucon if abs(tstat_polyjucon)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(pink%60) mlcolor(pink*1.5%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-20.1 70.1))   ///
	ylabel(-20(10)70, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)8, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "51") label(2 "Polyarchy > mean") label(5 "Jud. Constr. > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_poly_jucon.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_poly_jucon.eps", as(eps) replace


drop change_poly changese_poly tstat_poly

*** POLY & LEGCON
tab startyear_poly, gen(by_)
mrunning beta_poly exposure_poly flips_poly if !missing(beta_poly,beta_legcon, beta_polylegcon), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_poly
rename changese__1 changese_poly
drop change__* changese__*
gen tstat_poly = change_poly/changese_poly
drop by_*

tab startyear_legcon, gen(by_)
mrunning beta_legcon exposure_legcon flips_legcon if !missing(beta_poly,beta_legcon, beta_polylegcon), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_legcon
rename changese__1 changese_legcon
drop change__* changese__*
gen tstat_legcon = change_legcon/changese_legcon
drop by_*

tab startyear_polylegcon, gen(by_)
mrunning beta_polylegcon exposure_polylegcon flips_polylegcon if !missing(beta_poly,beta_legcon, beta_polylegcon), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_polylegcon
rename changese__1 changese_polylegcon
drop change__* changese__*
gen tstat_polylegcon = change_polylegcon/changese_polylegcon
drop by_*
* 64 countries 

twoway ///
	(hist exposure_polylegcon, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_poly exposure_poly, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(line change_legcon exposure_legcon, sort ///
		lw(.8) lcolor(gs7*1) lpattern(dash)) ///
	(scatter change_legcon exposure_legcon if abs(tstat_legcon)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_legcon exposure_legcon if abs(tstat_legcon)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(line change_polylegcon exposure_polylegcon, sort ///
		lw(.8) lcolor(emerald*1) lpattern(solid)) ///
	(scatter change_polylegcon exposure_polylegcon if abs(tstat_polylegcon)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(emerald*1.5%60) mlwidth(thin)) ///
	(scatter change_polylegcon exposure_polylegcon if abs(tstat_polylegcon)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(emerald%60) mlcolor(emerald*1.5%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-20.1 80.1))   ///
	ylabel(-20(10)80, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)12, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "64") label(2 "Polyarchy > mean") label(5 "Leg. Constraints > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_poly_legcon.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_poly_legcon.eps", as(eps) replace



***
*** Graph: Comparison interaction with liberal component
***

use "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_lib_0_diagnostics_4.dta", clear 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_lib
}
drop nwbcode					  
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_0_assoc_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_libassoc
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_0_exp_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_libexp
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_0_fairel_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_libfairel
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_exp_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_exp
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_assoc_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_assoc
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_fairel_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_fairel
}
drop nwbcode _merge

*** No horserace

*** LIB & ASSOC
tab startyear_lib, gen(by_)
mrunning beta_lib exposure_lib flips_lib if !missing(beta_lib,beta_assoc, beta_libassoc), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_lib
rename changese__1 changese_lib
drop change__* changese__*
gen tstat_lib = change_lib/changese_lib
drop by_*

tab startyear_assoc, gen(by_)
mrunning beta_assoc exposure_assoc flips_assoc if !missing(beta_lib,beta_assoc, beta_libassoc), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_assoc
rename changese__1 changese_assoc
drop change__* changese__*
gen tstat_assoc = change_assoc/changese_assoc
drop by_*

tab startyear_libassoc, gen(by_)
mrunning beta_libassoc exposure_libassoc flips_libassoc if !missing(beta_lib,beta_assoc, beta_libassoc), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_libassoc
rename changese__1 changese_libassoc
drop change__* changese__*
gen tstat_libassoc = change_libassoc/changese_libassoc
drop by_*

* 65 countries 

twoway ///
	(hist exposure_libassoc, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_lib exposure_lib, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs0*1%60) mlwidth(thin)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs0%60) mlcolor(gs0*1%60) mlwidth(thin)) ///
	(line change_assoc exposure_assoc, sort ///
		lw(.8) lcolor(gs0*1) lpattern(dash)) ///
	(scatter change_assoc exposure_assoc if abs(tstat_assoc)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_assoc exposure_assoc if abs(tstat_assoc)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs0%60) mlcolor(gs0*1.5%60) mlwidth(thin)) ///
	(line change_libassoc exposure_libassoc, sort ///
		lw(.8) lcolor(purple*1) lpattern(solid)) ///
	(scatter change_libassoc exposure_libassoc if abs(tstat_libassoc)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(purple*1.25%60) mlwidth(thin)) ///
	(scatter change_libassoc exposure_libassoc if abs(tstat_libassoc)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(purple%60) mlcolor(purple*1.25%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-5.1 50.1))   ///
	ylabel(0(10)50, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)14, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "63") label(2 "Liberal Comp. > mean") label(5 "Association > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_lib_assoc.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_lib_assoc.eps", as(eps) replace


drop change_lib changese_lib tstat_lib

*** LIB & EXP
tab startyear_lib, gen(by_)
mrunning beta_lib exposure_lib flips_lib if !missing(beta_lib,beta_exp, beta_libexp), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_lib
rename changese__1 changese_lib
drop change__* changese__*
gen tstat_lib = change_lib/changese_lib
drop by_*

tab startyear_exp, gen(by_)
mrunning beta_exp exposure_exp flips_exp if !missing(beta_lib,beta_exp, beta_libexp), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_exp
rename changese__1 changese_exp
drop change__* changese__*
gen tstat_exp = change_exp/changese_exp
drop by_*

tab startyear_libexp, gen(by_)
mrunning beta_libexp exposure_libexp flips_libexp if !missing(beta_lib,beta_exp, beta_libexp), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_libexp
rename changese__1 changese_libexp
drop change__* changese__*
gen tstat_libexp = change_libexp/changese_libexp
drop by_*
* 67 countries 

twoway ///
	(hist exposure_libexp, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_lib exposure_lib, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs0*1%60) mlwidth(thin)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs0%60) mlcolor(gs0*1%60) mlwidth(thin)) ///
	(line change_exp exposure_exp, sort ///
		lw(.8) lcolor(gs7*1) lpattern(dash)) ///
	(scatter change_exp exposure_exp if abs(tstat_exp)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_exp exposure_exp if abs(tstat_exp)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(line change_libexp exposure_libexp, sort ///
		lw(.8) lcolor(gold*1) lpattern(solid)) ///
	(scatter change_libexp exposure_libexp if abs(tstat_libexp)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(gold*1.5%60) mlwidth(thin)) ///
	(scatter change_libexp exposure_libexp if abs(tstat_libexp)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(gold%60) mlcolor(gold*1.5%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-15.1 40.1))   ///
	ylabel(-10(10)40, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)14, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "67") label(2 "Liberal Comp. > mean") label(5 "Expression > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_lib_exp.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_lib_exp.eps", as(eps) replace

drop change_lib changese_lib tstat_lib

*** LIB & FAIREL
tab startyear_lib, gen(by_)
mrunning beta_lib exposure_lib flips_lib if !missing(beta_lib, beta_fairel, beta_libfairel), gen(change__) ci gense(changese__)  adjust(by_*) nograph
rename change__1 change_lib
rename changese__1 changese_lib
drop change__* changese__*
gen tstat_lib = change_lib/changese_lib
drop by_*

tab startyear_fairel, gen(by_)
mrunning beta_fairel exposure_fairel flips_fairel if !missing(beta_lib, beta_fairel, beta_libfairel), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_fairel
rename changese__1 changese_fairel
drop change__* changese__*
gen tstat_fairel = change_fairel/changese_fairel
drop by_*

tab startyear_libfairel, gen(by_)
mrunning beta_libfairel exposure_libfairel flips_libfairel if !missing(beta_lib, beta_fairel, beta_libfairel), gen(change__) ci gense(changese__) adjust(by_*) nograph
rename change__1 change_libfairel
rename changese__1 changese_libfairel
drop change__* changese__*
gen tstat_libfairel = change_libfairel/changese_libfairel
drop by_*
* 67 countries 


twoway ///
	(hist exposure_libfairel, yaxis(2) w(5) start(0) freq fcolor(gs12%50) lcolor(white) barwidth(4.5)) ///
	(line change_lib exposure_lib, sort ///
		lw(.8) lcolor(gs0*1) lpattern(shortdash)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1%60) mlwidth(thin)) ///
	(line change_fairel exposure_fairel, sort ///
		lw(.8) lcolor(gs7*1) lpattern(dash)) ///
	(scatter change_fairel exposure_fairel if abs(tstat_fairel)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(white%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(scatter change_fairel exposure_fairel if abs(tstat_fairel)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(none)  msize(medlarge) mfcolor(gs7%60) mlcolor(gs7*1.5%60) mlwidth(thin)) ///
	(line change_libfairel exposure_libfairel, sort ///
		lw(.8) lcolor(orange*1) lpattern(solid)) ///
	(scatter change_libfairel exposure_libfairel if abs(tstat_libfairel)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(orange*1.5%60) mlwidth(thin)) ///
	(scatter change_libfairel exposure_libfairel if abs(tstat_libfairel)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(orange%60) mlcolor(orange*1.5%60) mlwidth(thin)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	ytitle("Country Count", size(small) axis(2))  ///
	yscale(range(-15.1 30.1))   ///
	ylabel(-10(10)30, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)12, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(1) order(2 5 8 - "Countries:" 1) ///
	label(1 "67") label(2 "Liberal Comp. > mean") label(5 "Fair Elections > mean") label(8 "Interaction")	 ///
	size(small)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_low_levels_st_interaction_lib_fairel.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_low_levels_st_interaction_lib_fairel.eps", as(eps) replace

exit 

***
*** Graph: Comparison interaction_0 with other cutoffs 
***

* ROLAW
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_neg25_rolaw_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_`kk'_rolaw_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(midblue*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(midblue*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(midblue*1) mlcolor(midblue*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-20.1 40.1))  ///
	ylabel(-20(10)40, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-20(10)40, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-20.1 40.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_poly_rolaw.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_poly_rolaw.eps", as(eps) replace


* JUCON
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_neg25_jucon_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_1949_2018_poly_`kk'_jucon_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(pink*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(pink*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(pink*1) mlcolor(pink*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-20.1 71.1))  ///
	ylabel(-20(10)70, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-20(10)70, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)70, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-20.1 71.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_poly_jucon.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_poly_jucon.eps", as(eps) replace



* LEGCON
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_neg25_legcon_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_poly_`kk'_legcon_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(emerald*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(emerald*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(emerald*1) mlcolor(emerald*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-20.1 80.1))  ///
	ylabel(-20(10)80, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-20(10)80, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-20.1 50.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_poly_legcon.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_poly_legcon.eps", as(eps) replace



* EXP
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_neg25_exp_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_`kk'_exp_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(gold*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(gold*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(gold*1) mlcolor(gold*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-20.1 50.1))  ///
	ylabel(-20(10)50, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-20(10)50, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-20.1 50.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_lib_exp.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_lib_exp.eps", as(eps) replace



* ASSOC
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_neg25_assoc_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_`kk'_assoc_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(purple*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(purple*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(purple*1) mlcolor(purple*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-10.1 73.1))  ///
	ylabel(-10(10)70, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-10(10)70, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-10.1 73.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_lib_assoc.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_lib_assoc.eps", as(eps) replace


* FAIREL
use "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_neg25_fairel_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_st_Inter_simple_1959_2018_lib_`kk'_fairel_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort ///
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///
		lcolor(orange*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(orange*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(orange*1) mlcolor(orange*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(6) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Democracy Effect (in %)", size(small))  ///
	yscale(range(-10.1 30.1))  ///
	ylabel(-10(10)30, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-10(10)30, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(small) axis(2))  ///
	yscale(range(-10.1 30.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(small) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_lib_fairel.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_lib_fairel.eps", as(eps) replace

